home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / QFLOAT / QSQRTA.C < prev    next >
C/C++ Source or Header  |  1996-03-24  |  2KB  |  107 lines
  1. /*    qsqrta.c        */
  2. /* Square root check routine, done by long division. */
  3. /* Copyright (C) 1984-1988 by Stephen L. Moshier. */
  4.  
  5.  
  6. #include "qhead.h"
  7. #include "mconf.h"
  8.  
  9. extern QELT qzero[], qone[];
  10. static QELT rndbit[NQ+1] = {0};
  11. static QELT rndset = 0;
  12. int mtherr(), qmovz(), shdn1(), shup1(), addm(), cmpm(), subm(), mdnorm();
  13.  
  14. int qsqrt( x, y )
  15. QELT *x, *y;
  16. {
  17. QELT temp[NQ+1], num[NQ+1], sq[NQ+1], xx[NQ+1];
  18. int i, j, k, n;
  19. long m, exp;
  20.  
  21. if( rndset == 0 )
  22.     {
  23.     qclear( rndbit );
  24.     rndbit[NQ] = 0;
  25.     rndbit[NQ-1] = 1;
  26.     rndset = 1;
  27.     }
  28. /* Check for arg <= 0 */
  29. /*i = qcmp( x, qzero );*/
  30. if( x[1] == 0 )
  31.     goto iszero;
  32. if( x[0] != 0 )
  33.     {
  34.     mtherr( "qsqrt", DOMAIN );
  35. iszero:
  36.     qclear(y);
  37.     return 0;
  38.     }
  39.  
  40. /* Bring in the arg and renormalized if it is denormal. */
  41. qmovz( x, xx );
  42. m = xx[1]; /* local long word exponent */
  43. /*
  44. if( m == 0 )
  45.     m -= normlz( xx );
  46. */
  47. /* Divide exponent by 2 */
  48. m -= EXPONE - 1;
  49. exp = (QELT) ( (m / 2) + (EXPONE - 1) );
  50.  
  51. /* Adjust if exponent odd */
  52. if( (m & 1) != 0 )
  53.     {
  54.     if( m > 0 )
  55.         exp += 1;
  56.     shdn1( xx );
  57.     }
  58.  
  59. qclear( sq );
  60. sq[NQ] = 0;
  61. qclear( num );
  62. num[NQ] = 0;
  63. n = WORDSIZE / 2;
  64.  
  65. for( j=0; j<2*NQ-6; j++ )
  66.     {
  67. /* bring in next word of arg */
  68.     if( j <= NQ-3 )
  69.         num[NQ] = xx[j+3];
  70. /* Do additional bit on last outer loop, for roundoff. */
  71.     if( j == 2*NQ-7 )
  72.         n = n = WORDSIZE / 2  + 1;
  73.     for( i=0; i<n; i++ )
  74.         {
  75. /* Next 2 bits of arg */
  76.         shup1( num );
  77.         shup1( num );
  78. /* Shift up answer */
  79.         shup1( sq );
  80. /* Make trial divisor */
  81.         for( k=0; k<=NQ; k++ )
  82.             temp[k] = sq[k];
  83.         shup1( temp );
  84.         addm( rndbit, temp );
  85. /* Subtract and insert answer bit if it goes in */
  86.         if( cmpm( temp, num ) <= 0 )
  87.             {
  88.             subm( temp, num );
  89.             sq[NQ-1] |= 1;
  90.             }
  91.         }
  92.     }
  93.  
  94. exp -= 1; /* Adjust for extra, roundoff loop done. */
  95. sq[1] = exp;
  96.  
  97. /* Sticky bit = 1 if the remainder is nonzero. */
  98. k = 0;
  99. for( i=3; i<=NQ; i++ )
  100.     k |= num[i];
  101.  
  102. /* Renormalize and round off. */
  103. mdnorm( sq, k );
  104. qmov( sq, y );
  105. return 0;
  106. }
  107.